arXiv:1502.01553vl [math.NA] 5 Feb 2015 


MINIMAL DEGREE il(curl) AND H{div) CONFORMING FINITE 
ELEMENTS ON POLYTOPAL MESHES 

WENBIN CHEN AND YANQIU WANG 


Abstract. We construct ff(curl) and ff(div) conforming finite elements on 
convex polygons and polyhedra with minimal possible degrees of freedom, i.e., 
the number of degrees of freedom is equal to the number of edges or faces of 
the polygon/polyhedron. The construction is based on generalized barycentric 
coordinates and the Whitney forms. In 3D, it currently requires the faces of the 
polyhedron be either triangles or parallelograms. Formula for computing basis 
functions are given. The finite elements satisfy discrete de Rham sequences 
in analogy to the well-known ones on simplices. Moreover, they reproduce 
existing H{cur\)-H(d'lv) elements on simplices, parallelograms, parallelepipeds, 
pyramids and triangular prisms. Approximation property of the constructed 
elements is also analyzed, by showing that the lowest-order simplicial Nelelec- 
Raviart-Thomas elements are subsets of the constructed elements on arbitrary 
polygons and certain polyhedra. 


1. Introduction 

On a contractible smooth manifold T C it is well-known Him [5] that the 
extended L? de Rham complex 

(1.1) 0-^ HAP{T) HA\T) -^ HK^{T) -^ 0, 

is exact, where d is the exterior derivative, and HK^{T), k = 0,... ,m, are Hilbert 
spaces containing all differential fc-forms w, such that both w and dcu are in L^. 
Using traditional vector proxy notation of differential forms, the de Rham complex 
can be expressed in 3D as 


0-^ H^{T) iJ(curl, T) R(div, T) 

• L^iT) 

-^0, 

and in 2D as either one of the following 




0-^ K H\T) H(curl, T) 

^ L\T) - 

-^0, 


0-^ K H\T) £r(div, T) 

^ L\T) - 

-^0, 



-a, 


d:r. 


where we conveniently denote the 2D curl operator by curl = 
the two complexes in 2D are indeed equivalent under the following mapping 

0 -1 


Note that 


il(curl, T) —>. R(div, T), where x = 


1 0 
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Thus it suffices to only study one of them, and in this paper we pick the one 
containing 77(div, T). 

The idea of finite element exterior calculus is to build finite dimensional sub¬ 
complexes of (O, and then patch the local discrete spaces on each mesh element, 
usually a polytope, together to obtain the finite element space on the entire mesh. 
To build conforming finite element spaces, certain continuity conditions will be 
imposed on the boundary of T. When T is a simplex or a hypercube, it is well- 
known that such sub-complexes can be built using polynomials, i.e., , V~A^ 

and T-LrA^, for 0 < /c < m (see [5] for definition of these spaces). Here we are 
interested in more general polygonal/polyhedral domain T, on which polynomial 
spaces like VrA^, V~ A^ and 'HrA^ are usually not enough for building conforming 
finite elements. For example, in 2D, one can not build TJ^-conforming, piecewise 
linear/bilinear, scalar hnite element space on meshes containing n-gons with n > 4. 
A solution is to use the generalized barycentric coordinates: Wachspress, Sibson, 
harmonic, and mean value, etc. (see [H na [201 m IMI IMl 113 Ea HE] and references 
therein), which allows one to build -conforming scalar finite element spaces using 
a larger set of basis functions [111 123 133 133 113 Ell ES E3 El]- For example, 
the Wachspress element uses rational functions. We would also like to mention 
two methods related to the generalized barycentric coordinates: the mimetic finite 
difference method (see the recent survey paper |32j ) and the virtual element method 
[44j . Both methods are defined on general polytopes. Among them, the lowest order 
virtual element method is indeed equivalent to an conforming finite element 
using a set of harmonic barycentric coordinates. 

Recall the traditional polynomial-valued barycentric coordinates defined on sim- 
plices, generalized barycentric coordinates {A^}, for i from 1 to the number of 
vertices, can be viewed as extensions of traditional barycentric coordinates to a 
polytope T. According to the construction, they may have some nice properties, 
which will be further explained later. In general, we expect {A^} to form a basis for 
an conforming scalar finite element on T. Extending such elements to iF(curl) 
and iJ(div) on general polytopes is not easy. As early as in 1988, researchers have 
realized the important role of Whitney forms in constructing vector-valued finite 
element spaces [Tj- The Whitney 1-form and Whitney 2-form on simplices are 
defined, respectively, by 

(1.2) = AVA, - A,VA, 

(1.3) Xj X VAfc -j- AjVAfc x VA^ -1- A^VA^ x VAj. 

Formally, by using generalized barycentric coordinates, they can be extended to 
general polytopes. There were several pioneering works on extending the Whitney 
forms and building iF(curl)/iJ(div) conforming finite elements over non-simplicial 
polytopes, including polygons [14], rectangular grids [23, and pyramids [23- In 
recent years, this idea has attracted more attentions. Gillette and Bajaj [23122] 
constructed dual mixed finite elements on polytopal meshs generated by taking 
the dual of simplicial meshes. Later in [3, Bossavit constructed edge-based and 
face-based Whitney forms on tetrahedra, hexahedra, triangular prisms, and pyra¬ 
mids using techniques called ‘conation’ and ‘extrusion’. And in the most recent 
work [23], Gillette, Rand and Bajaj constructed FF(curl) and iF(div) conforming 
finite elements on arbitrary polytopes using the span of all Whitney 1-forms and 
2-forms, respectively. We would also like to mention a few related works not using 


H{cur\) AND H(div) ELEMENTS ON POLYTOPAL MESHES 


3 


the Whitney forms. Kuznetsov and Repin [SOllST] constructed H{div) elements on 
polytopes with simplicial refinements by solving a local discrete mixed problem. 
Christiansen m constructed iJ(curl) and 7J(div) conforming finite elements on 
polytopes by using harmonic basis functions, which are known to be almost non- 
computable. Klausen, Rasmussen and Stephansen [29] directly constructed i/(div) 
conforming elements on polygons and simple polyhedra using generalized barycen- 
tric coordinates. A polyhedron in 3D is simple if all its vertices are connected to 
exactly 3 edges. The elements constructed in |2S]) although having minimal degrees 
of freedom, does not fit easily into a de Rham sequence. 

The main purpose of this paper is to provide a unified, easy-to-compute, and 
minimal degree construction of i?(curl) and H(div) conforming finite elements on 
convex polytopes, that satisfy the discrete de Rham sequence. Let us briefly explain 
how our work will be different from the existing results mentioned above. We aim 


at building sub-complexes of (1.11 using the minimal amount of basis functions that 


ensures iL(curl) and 7L(div) conformity. At the same time, we want the element 
to be constructed provides at least 0{h) approximation rate. Let us first recall the 
spaces constructed in [24]. Define 

WA°(T) = span{Xi}, WA\T) = span{W,j}, >VA^(T) = span{Wijk}- 

In [21], the authors have proved that the above defined finite element spaces are 
/H{curl)/H(div) conforming and contain h^{T), the lowest-order Nedelec- 
Raviart-Thomas spaces on simplices defined as following: 

In 2D: WA°(r) AV^AP{T) = span{l,x,y}, 

x(yV’A^(r)) D A^(r) = {ax -l- c, for a e R, c G K^}, 

(1.4) In 3D: WA°(r) DiPf A°(r) = span{l,x,y,z}, 

WAi(r) AiPf Ai(r) = (a X x-kb, for a,b G M.^}, 
WA^{T) 2VrA^{T) = (ax-kc, for a G K, c G K^}. 


Moreover, if T is a simplex, then WA^{T) coincides with Vi A^(T), i.e., all D in 
the above become =. 

Clearly, yVA°(T) is one of the smallest possible scalar finite elements on T that 
can ensure conformity. However, >VA^(T)/>VA^(T) are far from the smallest 
iL(curl)/iL(div) conforming elements on general polytopes. Indeed, denote by n 
the total number of vertices in T, then one has 


total number of Wij 
total number of Wijk 

For example, when T is a 3D cube, the above two numbers are 28 and 56, re¬ 
spectively. It is not clear whether Wij (or Wijk) are linearly independent or not. 
Thus one may need to use the least squares method in the implementation. Com¬ 
paring to the known smallest vector-valued finite element complex on a cube |36j . 
which uses 12 basis functions in the H(curl) element and 6 basis functions in the 
iL(div) element, the spaces WA^{T) and WA?{T) may contain too much redundant 
information. 
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We want to find the minimal discrete de Rham complex on general convex poly¬ 
topes that provides conforming approximations in , 77(curl) and -ff(div). Be¬ 
cause of the nice property of Whitney forms mm, we limit our searching in subsets 
of Wh}^{T). That is, we shall construct finite elements AiA^{T) satisfying 

AdA°(T) = WA°(T) and MK^iT) CWh^(T) for/fc=l,2. 

Now let us look at the smallest possible dimension of A4A^{T), for k = 1,2, on 
convex polytopes. We start from the 3D case. Denote by #R, and #F the 
number of vertices, edges and faces of a convex polyhedron T. Then, one has 
dimMA^{T) = dimWA°{T) = To ensure i/(curl) and H{div) conformity, 

which in turn requires tangential components and normal components be continuous 
across interfaces, respectively, our conjecture is that 

min (dimAdA^(T)) = min (dim AdA^(T)) = 

which remains to be verified later by construction. According to Euler’s formula 
for convex polyhedra, one has 

#F = #E + #F - 2 = (#E - 1) + (#F - 1). 


This helps to formulate an exact sequence that we aim to build: 

MA\T) 


(1.5) 0 


c A^A*^(T') grad 
dim—^Y 


dim—^E 


curl A’^A^(T') div 
dim—^F 


0 . 


Analogously, when T is a 2D polygon, we aim at building an exact sequence 


( 1 . 6 ) 


MA°{T) ^ x{MA\T)) div,^^ 

dim=#V dim=#E=i^V 


In the rest of this paper, we shall focus on constructing x{AiA^{T)) in 2D, as well 
as AIA^(T) and AiA^{T) in 3D, that make sequences (1.5)-(1.6) exact, and more 
importantly, allows one to build H (curl) and H (div) conforming finite element 
spaces. 

The rest of the paper is organized as follows. We briefly introduce the definition 
and properties of the generalized barycentric coordinates in Section Assump¬ 
tions on the polytope T and the generalized barycentric coordinates will also be 
stated in this section. Then, in Section]^ we construct iJ(div) conforming element 
x{AiA^{T)) for arbitrary convex polygons in 2D, which satisfies (1.6). Our formula 
is different from, and easier to compute in practice than the 2D formula given in 
|14) . although the resulting basis functions may be identical. Moreover, when the 
polygon satisfy certain shape regularity conditions, we prove the optimal mixed 
finite element a priori error. Numerical results are presented too. In Section]^ we 
construct i/(curl) conforming element AiA^{T) and H{div) conforming element 
AiA^iT) in 3D, which satisfy (1.5). The current construction only works for poly¬ 
hedra whose faces are either triangles or parallelograms. Examples show that our 
construction, as one unified formula, reproduces existing minimal degree finite ele¬ 
ments on tetrahedra, rectangular boxes, pyramids, and triangular prisms. We also 
construct finite elements on a regular octahedron, which has never been done be¬ 
fore. Moreover, for certain type of polyhedra, we prove that A^{T) C MA^{T), 
for fc = 0,1, 2, which will ensure the approximation property of AiA^iT). 
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2. Generalized barycentric coordinates and assumptions 

Let T be a convex polygon or polyhedron with n vertices denoted by v^, for i = 
1,..., n. The generalized barycentric coordinates are functions Xi, for i = 1,..., n, 
that satisfy: 

(1) (Non-negativity) All Xi, for 1 < i < n, have non-negative value on T; 

(2) (Linear precision) For any linear function L(x) defined on T, one has 

n 

L(vi)Ai(x), for all x S T. 

i=l 

The linear precision property is indeed equivalent to the combination of the follow¬ 
ing two properties: for all x € T, 

n n 

(2.1) ^Ai(x) = l, ^A,(x)vi=x. 

i=l i=l 

Different types of generalized barycentric coordinates have been proposed in both 
2D and 3D. Reader’s may refer to [151 (13 HOI UHl IMl [Ml ESI ESI EE] and references 
therein for more details. When T is a simplex, all generalized barycentric coordi¬ 
nates are identical, and they are equal to the traditional barycentric coordinates 
on simplices, which span the space of all linear polynomials. 

The spaces Af A^(T) that we plan to construct in this paper will be based on 
generalized barycentric coordinates. In the construction, we do require certain 
properties from generalized barycentric coordinates, which will be listed below as an 
assumption. We will also explain that the following assumption is not unreasonable, 
since there exist generalized barycentric coordinates that satisfy all terms in the 
assumption. But here we choose to list them as assumptions instead of limiting our 
interest to specific coordinates, in order to provide a more general setting. 

Assumption 1: There exists a set of generalized barycentric coordinates on T 
satisfying the following: 

• (Lagrange property) For all 1 < i, j < n, one has Ai(v^) = Sij, where Sij is 
the Kronecker delta; 

• (Trace property) In 2D, each Xi is piecewise linear on dT. In 3D, each 
Xi degenerates into a 2D generalized barycentric coordinate satisfying As¬ 
sumption 1 on each face ofT. 

• (Smoothness) For all 1 < i < n, one has Xi G C^(T). 

Remark 2.1. Assumption 1 is not unreasonable. It has been proved in |18j that 
all 2D generalized barycentric coordinates on convex polygons satisfy the Lagrange 
property and the trace property. In 3D, the Wachspress coordinates [15] and the 
mean value coordinates m have been defined and studied. The Wachspress coordi¬ 
nates satisfy the Lagrange property and the trace property on all convex polytopes 
[49) . The mean value coordinates have been proved to satisfy the Lagrange prop¬ 
erty and the trace property on convex polytopes whose faces are all triangular m- 
Both the Wachspress and the mean value coordinates are known to be in C°° in 
the interior of T and have unique continuous extension to dT. 

In 3D, we will need to impose an additional assumption on the convex polyhedron 
T, which basically requires each face of T must be either a triangle or a parallel¬ 
ogram. To explain the reason for such a restrictive assumption, we first list some 
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special properties of 2D generalized barycentric coordinates on triangles and paral¬ 
lelograms. Denote by | • | the length/area/volume of an edge/polygon/polyhedron, 
depending on the context. 

Lemma 2.2. Consider a triangle T with vertices 1 < i < 3, ordered counter- 

clockwisely. Denote the barycentric coordinates by Xi, 1 < i < 3. Their gradients 
VAi are two-dimensional constant vectors. We have 

(2.2) det [VAi VAj] = 2\T\ ’ 

for (*,j)e{(l,2), (2,3), (3,1)}. 

Proof. Denote by Ci the edge opposite to vertex v^, for 1 < i < 3. Clearly, VAi is 
a constant vector orthogonal to e^, pointing from towards Vj, and with length 
^pj|. Denote by 6ij the internal angle of T formed by edges and Cj. Then we 
have 

det [VA* VAj] = |VAi| |VAj| sin(7r — 9ij) 

2\T\ _ 1 

WW~W\' 

This completes the proof of the lemma. □ 


\^l\ \^j\ . ^ 

——r —^ Sin = 

2|T| 2|T| ^ 


Lemma 2.3. Consider a parallelogram T with vertices v^, 1 < * < 4, ordered 
counter-clockwisely. Denote the Wachspress coordinates on T by Xi, 1 < i < 4. 
Their gradients VA^ are two-dimensional vectors. We have 

det [VAi VA 2 ] -I- det [VA 3 VA 4 ] = , 

(2.3) 

det [VA 2 VA 3 ] -|- det [VA 4 VAi] = . 


Proof. Without loss of generality, denote the vertices of T, in counter-clockwise 
order, by Vi : ( 0 , 0 ), V 2 : (hi, 0 ), V 3 : (hi -|-fch 2 , h 2 ), V 4 : (hh 2 , h 2 ), where hi, h 2 and 
k are positive constants. Then, one can easily compute the Wachspress coordinates 
and their gradients: 

,-{h2-y) X - 2 ky - hi-\-kh2 u 


Xi = 

X 2 = 
^3 = 
A4 = 


{hi-x + ky){h 2 - y) 

VAi 

hih2 


{x - ky){h 2 - y) 

VA 2 

hih 2 


{x - ky)y 

VA 3 

hih2 

{hi-x + ky)y 
hih2 

VA 4 


hih2 hih2 

1^2 —y —x-\-2ky — kh2^i 
^ hih 2 ’ h)h^ ^ ’ 
r y X - 2ky ^ 

^hiha’ hiha 

r -y 

^hiha’ hih 2 


The lemma hence follows from direct calculation. 


□ 


Now we state the additional assumption on T: 

Assumption 2: In 3D, assume each face of polyhedron T be either a triangle 
or a parallelogram. Moreover, assume the trace of the generalized barycentric coor¬ 
dinates chosen in our construction satisfy equations {2.2)-(2.3] on the faces ofT. 
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Remark 2.4. Equations (2.2|-(2.3) will later ensure that each function in the con¬ 
structed i/(div) finite element space has constant normal components on faces. 
This is why we need Assumption 2. Similar but much more complicated equations, 
with non-constant right-hand sides, can be obtained for general polygons. Whether 
they can be used to build vector-valued finite elements on polyhedra not satisfying 
Assumption 2 is a topic for future research. 


Remark 2.5. According to lemmas |2.2|2.3[ for convex polyhedra with only triangu¬ 
lar faces, both the Wachspress and the mean value coordinates can be used in the 
construction; while for convex polyhedra with both triangular faces and parallelo- 
gramal faces, only the Wachspress coordinates can be used. 


Throughout the rest of this paper, we always assume the polytope, as well as 
the generalized barycentric coordinates defined on it, satisfy Assumptions 1-2. It 
is known that all polygons and many polyhedra, including the most frequently 
used tetrahedra, parallelepipeds, triangular prisms, and pyramids, have generalized 
barycentric coordinates defined on them that satisfy these assumptions. 


3. Construction in 2D 

Let r be a convex polygon. Denote by v^, 1 < i < n, the vertices of T ordered 
counterclockwisely, and by the edge connecting vertices and v^+i, where we 
conveniently denote Vj = Vj i^„iod n) when the subscript j is not in the range of 
{1,..., n}. Similar tricks of indexing will be used frequently without special men¬ 
tioning. Denote by and the unit outward normal and the unit tangent vector 
in the counterclockwise orientation on . Choose an arbitrary point x» inside poly¬ 
gon T, and denote by R the triangle with base and apex x*. Denote by di the 
distance from x* to Ci. Let \ei\, \Ti\ and |T| be the length of e^, the area of R and 
T, respectively. It is clear that \R\ = and |T| = We use the stan¬ 

dard notation LP{T), W^’P{T), H^{T) and iL(div,T), with s S K and 1 < p < oo 
for different type of Sobolev spaces, equipped with corresponding innerproducts 
and norms. For simplicity, denote by || • Hr and || • ||e^ the norm on T and 
respectively, while by || • || i^t the norm on T. Finally, denote by hx the diameter 
of T. 


3.1. Discrete spa ce and basis function. Recall that RihP{T ) = sponjAi, i = 
I,..., n}. By (2.1), one has K C MA^{T) and thus the sequence ( |l. 6 [ ) is obviously 
exact at the MA^{T) node. In order to ensure the exactness at the x(AtA^(T)) 
node, we would like to define x(AtA^(T)) with an orthogonal decomposition, i.e., 
the discrete Helmholtz decomposition: 

x{MA\T)) = curlMA°(T) 0 (div^')®. 


where div^ stands for a pseudo-inverse of div under proper choice of spaces such 
that (div^)M contains functions orthogonal to curlAIA°(T) and with divergence in 
K. In practice, it is much easier if one relaxes the orthogonality a little bit through 
replacing 0 by A, and thus we consider the following construction: 

x{M.A^(T)) = curlAIA°(r) 0 span{x — x*} 

= spanjcurlAi, i = I,..., n} A span{x — x*}. 


Later we shall show that the above definition is independent of the choice of x*. 
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By construction, it is clear that div x(A^ A^(r)) = K and curlAdA°(r)nspan{x— 
X*} = {0}. Therefore, the sequence (1.6) is also exact at the x(AlA^(T)) node. 
Now we know that the entire sequence ( |1.6[ ) is exact. By counting dimensions and 
since obviously dimMhP{T) = n, one must have dimx{M.h}{T)) = n. Next, we 
explicitly construct a set of basis for x(AlA^(T)). 

For 1 < i, Z < n, define 


The above notation can be extended to indices not in {l,...,n} using modular 
arithmetic. 


Lemma 3.1. For each 1 <i <n, define G x(AlA^(T)) by 

71 

(3.2) qi = Cj,o(x - X*) + ^ Ci,fcCurlAfc, 

fc=i 

where Cj,o = ^ and d^k = - ^ h,k+i- Then, one has q^ • le^ = S^j for all 

1 ^ j < H., and the set {q^, 1 < f < n} form a basis for x(AlA^(T)). 


Proof. Notice that for all 1 < fc < n, one has 


n n " IT I 

Kk+i= X! ^ =Q’ 

1^1 1^1'' 

which implies that 

— (0 — ^i,k- 

n 

Therefore, by the definition of generalized barycentric coordinates and Assumption 
1 , we have 


^i,k Q,fc+1 — 



Qi ■ tijUj — Cifiijx. X*) • nj\ej + 'y ( Ci^feCurlAfe • njlcj — Ci,o dj ^ ( Ci^k q. 


= Ci,o- 

= 6ij. 




k=l 

led led 


fc=i 


dtn 


-Tn{]f\\TA + hj 


The set {q^, 1 < j < n} is linearly independent, because X]r=i “*1* ~ ^ implies 
that 0 = (X]”=i ai^i) ■ rijlej = aj for all 1 < j < n. Since dim x(AlA^(r)) = n, the 
set {q., 1 < i < n} must form a basis for x(AlA^(r)). This completes the proof of 
the lemma. □ 


Remark 3.2. From the basis it is clear that for any function q € x (7VfA^(T)), the 
normal component q • n is piecewise constant on dT. Moreover, the normal compo¬ 
nents on edges form a unisolvant set of degrees of freedom for x (A1A°(T)). Such 
a choice of degrees of freedom guarantees that one can build H(div) conforming 
finite element spaces on general polygonal meshes using x {MA^{T)). 
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Remark 3.3. When T is a triangle, all currently known generalized barycentric co¬ 
ordinates degenerate to the unique triangular barycentric coordinates {Ai, A 2 , A 3 }. 
In this case, the space A^A°(T) is identical to span{l,x,y}, and consequently the 
space X (A^A^(T)) is identical to Vi A^{T), the lowest-order Raviart-Thomas finite 
element on triangles. When T is a rectangle and A^’s are chosen to be the Wach- 
spress coordinates, the space AIA°(T) is identical to span{l,x,y,xy}, and con¬ 
sequently x(AdA^(T)) is identical to 'HiA^(T), the lowest-order Raviart-Thomas 
finite element on rectangles. In this sense, the space y (A^A^(T)) can be viewed 
as the extension of the lowest-order Raviart-Thomas finite element to general poly¬ 
gons. 

A more important relation between 7^j“A^(T) and y (A^A^(T)) is given in the 
following lemma: 


Lemma 3.4. 


The space y (A^A^(r)) reproduces all funetions in Vi A^{T), i.e., 
V^A\T)Cx{MA\T)). 


Proof. This follows immediately from the definition of y (A ^A^(r)) and the fact 
that VfA^{T) C Af A°(r), which comes from Equation ( |2.l[ ). □ 


Remark 3.5. Lemma 
the choice of x*. 


3.4 


indicates that the space y(AfA^(T)) is independent of 


Finally, we briefly show that y(Af A^(r)) constructed in this section is a subspace 
of y(yV’A^(T)). By Equation (2.1), it is not hard to see that 

n 

(3.3) ^ Wij = -VA„ for all 1 < z < n, 

1=1 

which implies that curlA4A°(r) = y(VAf A°(T)) C y(yV’A^(r)). Recall the inclu¬ 
sion relation of finite elements in (1.4), one has 7^]"A^(r) C y(>VA^(T)). Combining 
the above with the definition of y(Af A^(T) gives y(AfA^(r) C y(yV’A^(T)). 


3.2. Interpolation operator and its properties. To make sure that the mixed 
finite element theory works on the finite element y (At A^ (T)), we define an interpo¬ 
lation operator into y (AtA^(r)) which satisfies certain stability and approximation 
properties. For convenience, we introduce the notation <, > and « for ‘less than 
or equal to’, ‘greater than or equal to’, ‘both less than or equal to and greater than 
or equal to’ up to a constant independent of the shape of all polygons in a given 
mesh. 

Clearly, to establish any kind of stability and approximation properties, the 
polygonal mesh needs to satisfy certain shape regularity conditions. We assume for 
all polygons in the mesh, 

• The area of the polygon is related to its diameter as follows: 

\T \« 4; 


• The gradient of Ai, for 1 < i < n, on T satisfies 

(3.4) \yK\<hf\ atallxeT, 

where | • | stands for the Euclidean length. It has been proved in [19] that 

(3.4) holds for Wachspress coordinates as long as h*, the minimum distance 
from any vertex of T to a non-incidental edge, satisfies h* pe Ht- 






10 


WENBIN CHEN AND YANQIU WANG 


• The following trace inequality and approximation property of projection 
hold on T : 

ll</>llWT)</i?'ll</'llT + /*T||V</.|||, iovcj,€H\n 

U-PTnT<hT\mi,T, iovcj>&H\T), 

where Pt denotes the orthogonal projection onto M. It is known that 
when T satisfy certain shape regularity conditions, (3.51 holds on T. Read¬ 
ers may refer to uni EH uni on for further discussion. 

Now let use define the interpolation operator. For any q S H{div, T) n 
with p > 2, define Il-rq £ X (A^A^(T)) by 

(Erq) ■ I q • rij ds, for 1 < j < n. 

I Jej 

The requirement p > 2 is to guarantee that /g . q • nj ds be well-defined. One may 
circumvent this requirement by using Clement type interpolations [12) . According 
to the definition, it is clear that 

n ^ „ 

where at = —r / q • ds. 

Moreover, by the unisolvancy of the degrees of freedom and Lemma 3.4 we know 
that n-r preserves all functions in h.^{T), i.e.. 


rij" 


cx + a 
cx + b 


cx + a 
cx + b 


for all a, 6 , c £ 


Denote by It the nodal value interpolation into A4A°(T). Properties of nodal 
value interpolation for generalized barycentric coordinates have be discussed in 
[Mils]. Then we have: 


Lemma 3.6. Let p > 2. For any q £ iL(div, T) n (LP(T))^ , one has divIlTq = 
Prdivq. For any (p £ one has Hrcurlt/) = cmllTp- In other words, the 

following diagram is commutative: 


w^’P{T) p(div, T) n (LP(r))2 l^{t) 


It 

MA°iT) 


X{MK\T)) 



Proof. Given q £ iL(div, T) n {Lp{T)Y, let /g. q • ds for 1 < z < n. 

Then by the definition of basis function q^, one has 


n 

divH'rq = div 


E 




iii ds 


divqda; = Prdivq. 
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Given (j) G W^’^{T). Then /pt/) = Note that for all 1 < j < n, one 

has 


(n-pcurl^) • rijle^ = —/ curl;/) • ds 

Pi I Jei 


eil Je, dtj 


d(t> ^ - </>(vi+i) 


and 


d\i 


(curl/T0) • Ujle^ = ( curl^(^(v,)A, j • nj\e^ = - ^ 


2=1 




9t,- 


*'1 '^*'1 I'^ji 

By the unisolvancy of the degrees of freedom, we have IlT’Curlc^ = cut\It4’- This 
completes the proof of the lemma. □ 

To prove the stability and approximation properties of IIt’, we first derive the 
following estimate of q^: 

Lemma 3.7. For 1 < i < n, one has 

l|q*llT < G(n)|ei|, 

where C{n) is a general positive constant depending only on n. 

Proof. Note that 


1 ^ 


Iqillr = ||cj,o(x- X*) + y^Ci,fcCurlAfc||^ 
fc=i 

/ n \ 

< (n + 1) f g||x — x*|jy + ^ Ci,fc|j VAfclly 

71 

= (n + 1 )( Jo + Jfc). 


fc=i 


k=l 


For Jo, we have 


Jo = 




4|T|2 


|x - x^IIj. < 




m 


:\T\hl<\e.\\ 


Here in the last step we used the assumption |T| « h"^. Next, by (3.41, we have the 
following estimate for J^: 


Jk — Ci^fell VAfelly < c^.kj^ ^ Cj fc — ( ^ ^ 


1 


1=1 


< 


n —1 

2 i<r< 


[<l<n ' ' ' J 


Combining the above, we have proved the lemma. 


□ 


Denote by Qt the (L^(r))^ projection onto K^. Clearly we have WtQt^ = Qt<\- 
Next we prove the following technical lemma: 
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Lemma 3.8. For q € (i7^(T))^, one has 

||nT(q — QTq)||T ^ C'(?T-)^T||q||i,T, 
where C{n) is a general positive constant depending only on n. 

Proof. For convenience, denote q = q — Qt^I- Then by the Schwarz inequality and 
Lemma 13.71 


nTqllT = II XI ( / q-nids) q*||T <( Til / ^ ) llq*llr 


1 


^ llqlleJlqillT ^ Ill~ll2^ 

< (|e,|||qLJ . 


ed 


Then, by (3.5), one has 

llqlle. < h-^ 

Combining the above gives 


T llqllr ■ 


dT||Vq|l| < hrWqWl 


T- 


linTqllT c'N ( X l®*l ] ^ c'^dTlIqlli.-z 


\i=l / 

This completes the proof of the lemma. 

Next we prove the following stability property of 11^: 
Lemma 3.9. For q € (H^(T))^, one has 

l|nTq|l//(div,T) ^ C'(n)||q||i,r- 


□ 


Proof. By Lemma 3.6 it is clear that we o nly need to prove Ilfl-rqllT C'(n)||q||i_r. 
Using the triangle inequality, Lemma 3.8 and the stability of the Lf projection Qt, 
one has 

lin^qllT < linT(q — QTq)l|T + linT<3Tq|lT 

< C'(n)/iT||q||i,T + IIQrqllT 

< ^(njllqlli^T. 

In the above we have used the fact that = QtH- This completes the proof 

of the lemma. □ 

Finnally, we prove the approximation property of IIt’: 

Lemma 3.10. For all q S {Fl^{T)Y, one has 

||q — n-rqllr ^ C'(n)/iT|lq|li,T- 
Moreover, z/divq S F[^{T), then one has 

||div(q- riTqjllr < /iT||divq|| i,t. 

Proof. By the triangle inequality, the fact that nTQrq = Qrq, Lemma [TSl and 
the approximation property of Qt (similar to (3.51), one has 

||q — Ilrqllr (5; ||q — Qrqllr + ||nT(q — QTq)||r 
< C{n)hT\\(l\\i,T- 


The second part of the lemma follows from Lemma 3.6 and Inequality (3.5). □ 
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Remark 3.11. Because of the above properties of Hr, the finite element x (A1(T)) 
fits the theoretical framework of mixed finite element methods in the book by 
Brezzi and Fortin |^, as long as the polygonal mesh satisfies all shape regularity 
assumptions and the number of vertices in each polygon is bounded above. In 
this case, the mixed finite element achieves optimal approximation error in both 
II ■ l|L 2 (n) and ||div(-)||i 2 (f 2 ), where denotes the entire computational domain. 


3.3. Numerical results. In this section, we first draw a set of basis {q^} for 
i?(div) element on a random pentagon in Figure]^ in order to give the reader a 
direct picture of these basis functions. The basis is generated using the formula 
(3.2), with Xi set as the Wachspress coordinates. 



Figure 1. Basis {q^} for iJ(div) element on a random pentagon. 



Figure 2. Meshes of size 8x8. (1) A quadrilateral mesh. (2) A 
hexagonal mesh, with mostly hexagons and a few pentagons and 
quadrilaterals. It is generated as the dual mesh of an 8 x 8 uniform 
triangular mesh, as shown in dotted lines. (3) Centroidal Voronoi 
tessellation consisting of 8 x 8 cells (see m and references therein). 

Consider the Poisson’s equation on (0,1) x (0,1) with Dirichlet boundary con¬ 
dition. We test this problem on three different types of meshes, as shown in Figure 
Wachspress coordinates are used to define Xi. The example problem is solved on 
a sequence of meshes, using the mixed finite element method with x(Af A^(T))-M 
discretization. Denote by p and u the exact flux and the exact primal solution, 
while by Ph and Uh the corresponding numerical solutions. We first set the exact 
solution to be M = sin(7rx) sin(7ri/), which is smooth. Numerical results are reported 
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in Tables [T|^ in which the ‘order’ is the value of r in 0{h^) computed using the 
errors on two consecutive meshes. From the table we can see that ||p — Ph\\L^, 
||divp — divp;i||i 2 and ||u — Uh\\L^ have at least 0{h) convergence, which agrees 
well with the theoretical prediction. We also point out that although the centroidal 
Voronoi tessellation in Figurej^appears to contain very short edges, which may the¬ 
oretically break the condition given in |19| for the assumption (3.4), the numerical 
results presented in Table seem to be unaffected. 


Table 1. Example problem with exact solution u = 
sin(7ra;) sin(7r2/). Errors of the x(WlA^(T))-]R approximation 
on quadrilateral meshes as shown in Figure 


Mesh Size 

Up - Phh^ 

error order 

lldivp - divp/i||i 2 
error order 

||W - M/i||l2 

error order 

4x4 

5.2843e-l 

3.1580e-b0 

1.6184e-l 

8 x8 

2.6040e-l 

1.0210 

1.6087e-b0 

0.9731 

8.1764e-2 

0.9850 

16 X 16 

1.2971e-l 

1.0054 

8.0813e-l 

0.9932 

4.0974e-2 

0.9968 

32 X 32 

6.4810e-2 

1.0010 

4.0454e-l 

0.9983 

2.0498e-2 

0.9992 

64 X 64 

3.2405e-2 

1.0000 

2.0233e-l 

0.9996 

1.0251e-2 

0.9997 

128 X 128 

1.6204e-2 

0.9999 

1.0117e-l 

0.9999 

5.1255e-3 

1.0000 

256 X 256 

8.1023e-3 

0.9999 

5.0587e-2 

0.9999 

2.5628e-3 

1.0000 

512 X 512 

4.0513e-3 

0.9999 

2.5293e-2 

1.0000 

1.2814e-3 

1.0000 


Table 2. Example problem with exact solution u = 
sin(7ra;) sin(7ry). Errors of the x(WlA^(T))-K approximation 
on hexagonal meshes as shown in Figure 


Mesh Size 

Up - Phh^ 

error order 

lldivp - divp/i||i 2 
error order 

||w - ■U/i||l2 
error order 

4x4 

2.7502e-l 

2.6008e-b0 

1.3488e-l 

8 x8 

1.0994e-l 

1.3228 

1.4988e-b0 

0.7951 

7.6665e-2 

0.8150 

16 X 16 

4.5041e-2 

1.2874 

7.9379e-l 

0.9170 

4.0330e-2 

0.9267 

32 X 32 

2.0013e-2 

1.1703 

4.0721e-l 

0.9630 

2.0646e-2 

0.9660 

64 X 64 

9.4150e-3 

1.0879 

2.0608e-l 

0.9826 

1.0442e-2 

0.9835 

128 X 128 

4.5673e-3 

1.0436 

1.0365e-l 

0.9915 

5.2510e-3 

0.9917 

256 X 256 

2.2498e-3 

1.0215 

5.1973e-2 

0.9959 

2.6330e-3 

0.9959 

512 X 512 

1.1166e-3 

1.0107 

2.6023e-2 

0.9980 

1.3184e-3 

0.9979 


It would be interesting to compare the numerical results on quadrilateral meshes 
given in Table with the numerical results of the lowest order Raviart-Thomas 
element presented in [T]. The Raviart-Thomas element can be extended to con¬ 
vex quadrilaterals via the Piola transform associated to a bilinear isomorphism, 
but with a degeneration of approximation rate in ||div(p — p;t)||L 2 (see m)- It is 
not hard to check that, on quadrilaterals that are not parallelograms, the space 
X (AIA^(T)) is indeed different from the polynomial-valued, lowest-order Raviart- 
Thomas element via Piola transform, because in this case x(AIA^(T)) consists 
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Table 3. Example problem with exact solution u = 
sin(7ra;) sin(7ry). Errors of the approximation 

on centroidal Voronoi tessellations as shown in Figure 


Mesh Size 

Up - VhWh^ 

error order 

II divp - divp/i||i 2 
error order 

||w - M/i||l2 
error order 

4x4 

4.5335e-l 

3.1186e-H0 

1.6102e-l 

8 x8 

1.8368e-l 

1.3034 

1.5915e-h0 

0.9705 

8 .1220e-2 

0.9873 

16 X 16 

7.4684e-2 

1.2983 

7.7831e-l 

1.0320 

3.9513e-2 

1.0395 

32 X 32 

2.9515e-2 

1.3394 

3.9116e-l 

0.9926 

1.9829e-2 

0.9947 

64 X 64 

1.3361e-2 

1.1434 

1.9703e-l 

0.9893 

9.9831e-3 

0.9901 

128 X 128 

6.3094e-3 

1.0825 

9.7955e-2 

1.0082 

4.9627e-3 

1.0084 

256 X 256 

3.0048e-3 

1.0702 

4.8807e-2 

1.0050 

2.4726e-3 

1.0051 


of rational functions. Therefore, the x(A^A^(r))-K discretization will still pro¬ 
vide optimal 0{h) convergence rate in ||div(p — p;i)||L 2 , as shown in Table In 
comparison, numerical results given in [T], using the lowest order Raviart-Thomas 
element via Piola transform, does not convergence in ||divp — divp/j ||^2 when the 
mesh consists of general quadrilaterals. 


Table 4. Example problem with exact solution u € Errors 

of the A^(r))-K approximation on quadrilateral meshes as 
shown in Figure 


Mesh Size 

IIp-p/i||l2 
error order 

||w - UhWl-^ 
error order 

4x4 

9.1202e-2 

4.6653e-2 

8 x8 

6.5317e-2 

0.4816 

2.3850e-2 

0.9680 

16 X 16 

4.6480e-2 

0.4909 

1.2045e-2 

0.9856 

32 X 32 

3.2970e-2 

0.4955 

6.0512e-3 

0.9931 

64 X 64 

2.3350e-2 

0.4977 

3.0326e-3 

0.9967 

128 X 128 

1.6524e-2 

0.4989 

1.5181e-3 

0.9983 

256 X 256 

1.1689e-2 

0.4994 

7.5946e-4 

0.9992 

512 X 512 

8.2669e-3 

0.4997 

3.7984e-4 

0.9996 


We also test a second example problem, under the same settings but with exact 
solution u = ■sj\{p — x) — \p^, where p is the radius in polar coordinates. One 

can easily verify that —Au = 1 on (0,1) x (0,1), and moreover, u G 1)^). 

Numerical results for the second example problem using the quadrilateral meshes 
are reported in Table Note that ||divp — divp/jH^^a is not included since for 
this test problem, one has divp = divp^, = — 1. From the table, we observe that 
Up ~ is of approximately which is reasonable because p G 

while \\u — Uh\\L^ « 0{h) because u is in 
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4. Construction in 3D 


4.1. Definitions and properties. Let T be a convex polyhedron satisfying As¬ 
sumptions 1-2. Denote by v^, i = 1,... ,n, the vertices of T. Then, for each pair 
of indices {i,j}, 1 < *, j < n, we have the Whitney 1-form Wij. Similarly, for each 
triplet of indices {i,j,k}, 1 < i,j,k < n, we have the Whitney 2-form Wijk- It is 
not hard to see that Whitney forms have the following properties: 

= 0 , 

lAyfc = 0, if at least two of i,j, k are identical. 




Moreover, using the definition of Whitney forms, Equation (2.1) and elementary 
vector calculus identities, one has 


(4.1) curl = 2VA, x VA^ = 2 ^ Wijk- 

k=l 

We also state a result from [21]. Denote by = Vj — for all 1 < i, j < n. For 
any constant vector a S one has 

^ ^ (a • Tij)W,j = ^(a • Tij)W,j = a, 

i<j 

(4.2) ^ ^ ((a X Vi) • Tij)Wij = ^((a x v,) • Tij)Wij 

i<j 

= ^ = a X X. 

i<j 

The reason that Whitney forms are so important in the construction of H (curl) 
and iL(div) spaces is that, they naturally satisfy certain conditions on edges/faces 
of T. Before summarizing these in lemmas, we first need to clarify the concept of 
‘edges’. Denote by Cij the directed line segment pointing from Vi to Vj, and by \eij\ 
its length . Notice that Cij may not be a natural edge of polyhedron T. Indeed, we 
classify all Cij into three disjoint categories: 

(1) £ is the set of all Cij that coincides with a natural edge of T; 

(2) is the set of all lying on dT but not in £; 

(3) Si is the set of all in the interior of T, i.e., not lying on dT. 

An illustration of these categories is given in Figure [^ We point out that each 
category actually contains both and e^i, for a given pair of indices i and j. The 
union of all three categories covers all e^, for I < <n. Notice that Sp and Si 

can be empty for certain polyhedra. On each tij , denote by t^- the unit tangential 
vector pointing from v, to Vj. We emphasize that only the e^- G S will be called 
an ‘edge’ of T, while the others are just called ‘directed line segments’. 


Lemma 4.1. Let eu G £. Then for all I < *, j < n, one has 

( \^\ = ^kl , 

Wij ■ tki\ehi — \ ~\^\ */ 

I 0 otherwise. 


Proof. The proof follows immediately from the definitions of A^, Wij and Assump¬ 
tion 1, which states that A^ is linear on all Cki € S. □ 
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Figure 3. Illustration of three categories: ei2 € £, 652 € £f, es2 € £i- 


Remark 4.2. On eki G £f or £ 1 , we do not have results similar to Lemma |4.H since 
Xi may not even be linear on eu- 

Next we define another important form on each Cij G £. Denote by Fij the set of 
two faces of polyhedron T that share the edge eij , and by Vy the set of all vertices 
on Fij. For a fixed index 1 < i < n, note that any for eik G £i can be written 
as a linear combination of all for e^- G £. Such a linear combination is not 
uniquely defined if vertex is connected to more than 3 edges of the polyhedron. 
Nevertheless, we can always fix a linear combination for each vertex v^, and denote 
this chosen one by 

(4.3) ^ 

j,eijeS 

Now, define 

= w^j + li E - E 

\vkeVij,eikeSF VfeGVij. ejfcGfj? 

+ 5 1 E cSm.- E 

\k,eikeSi k,ejk&£i 

In the above, one may view IFy + 4 (Ev,gv,,. e,,G£i. “ Ev,GV.,.e,.G£p 

as the ‘surface’ component of IF^ and 4 (Efe.e,,G£z - J2k,e,kG£i 

as the ‘interior’ component of Wij. An illustration of the surface component of Wij, 
which can also be written as W^j + 4 (Ev,GV.,,e.fcG£p + Ev,GV.,.e,,G£i. ^kj) , 

is given in Figure Note that if both faces sharing are triangles, the surface 

component of IF^ is just Wtj. 

The vector function Wij has many nice properties. First, it is obvious that 
Wij = —Wji- Now, let us fix a direction for each edge of T. The collection of all 
edges in £, with the prefixed direction, is denoted by Similarly, one may denote 
the collection of all edges in £ with direction opposite to the prefixed one as £~. 
The two sets f + and £~ contain the same edges, but with opposite directions. For 
any two edges and Cki in £'^, denote by Seij^eki the Kronecker delta whose value 
is 1 if Cij = Cki and 0 otherwise. Then, we have the following lemmas: 

Lemma 4.3. The set {Wj, for e^ G £+} satisfy Wj ■ tki\eki = j^\Se,j,eki for all 
eki G £'^, and hence is linearly independent. 
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Figure 4. Illustration of surface component of Wij when e^- S 
8 is shared by two faces of T which are: (1) two triangles; (2) 
one triangle and one parallelogram; (3) two parallelograms. Here 
we conveniently use thick arrow to denote Wki and thin arrow to 
denote ^Wki on any eu- 


Proof. This follows immediately from the definition of Wij, Lemma 
fact that J2eijeS+ ^ implies that Cki = \eki\ {j2eiie£+ 

for all Cki & 8^. 


4.1 


and 

^kl\eki 


the 
= 0 
□ 


Lemma 4.4. It holds that K^lT) C span{Wij, for etj G 8^}. 


Proof. Let us first point out that span{Wij, for eg G 8'^} = span{Wij, for eg G 
8 }. By the definitions of Wg 
that V, X T 


and Cg, Equation (4.2), Assumption 2, and the fact 


for any a G one has 


((a X V,) • Tg)Wg = ((ax V,) • Tg)Wg 

Gij GS S-ij 

^ n 1 ^ 

((a X Vi) • Tik)Wik + 2 E E ((a X Vj) • Tjk)Wjk 

i—lk,eikGSF 3 — '^k,ejkGSF 

^ n 1 ^ 

+ 2 E E ((a X Vi) •'ri/c)Wifc + -^ ^ ((a X Vj) • Tjfc)IF,fc 

i—lk,eikGSi j—'^k,ejhGiSi 

n / ■JT' 

= X! (E((^ ^ ■ '^^k)Wik 

i=i \fe=i 

= 2a X X. 


This indicates that a x x G span{Wg, for eg G 8^}. Similarly, one can prove that 
for any b G K^, 

(b • Tg)Wg = 2b. 

Recall that Vi h.^(T) = span{sL x x + b, for all a, b G M^}. This completes the 
proof of the lemma. □ 


Denote by F the set of all faces of T, and by n/ the unit outward normal vector 
on f G F with respect to T. For each f G F, denote by |/| its area and by df the 
oriented boundary of / such that its orientation satisfies the right-hand rule with 
Uf. If eg lies on df and has the same direction as the orientation of df, we say 
eg G df. If eg lies on df and has the opposite direction as the orientation of df, 
we say eg G —df. 
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Lemma 4.5. Let f € and Cij € £, then one has 


curlWi- 


n/l/=<!-^ 


l)y if Cij e 5/, 

if Cij G -df, 
0 otherwise. 


Proof. Notice that for any f G iF and 1 < i,j < n, by Equation (4.1), one has 
cuTlWij ■ nf\f = 2(VAi X VAj) • nf\f = 2(V/Ai x V/Aj) • nf\f, 

where VyAi|/ denotes the tangential component of VA^ on /. By Assumption 1, 
VfXi\f is non-zero only if is a vertex on face /. It is then clear that curlVEy • 
nf\f is non-zero only when both and Vj are vertices of face /. Consequently, 
cmlWij ■ nf\f is non-zero only when Cij G df or —df. 

For f G P, denote by V(/) the set of vertices on face /. Without loss of 
generality, assume / lies on the xy-plane with outward normal rif = [0,0,1]‘, and 
denote by A^^\ for all v^, G V(/), the 2-dimensional barycentric coordinates on 
polygon /. By Assumption 1, the 3D coordinate Afc, where G V(/), degenerates 


( 2 ) 

to A^ ' on /. Consequently, V/Afc is equal to 


V(2)Af)' 
0 


, where stands for the 


2-dimensional gradient on the cry-plane. Note we have for all eu G df that 


(V/A/c X V/A;) • nf\f = 


( 

■v(2)a^2)' 

X 

■V(2);)^(2)' 

j. 

O' 

0 


0 


0 

/ 

1 




Consider the case Cij G df. When / is a triangle, it is clear that by Lemma 2.2 


curlWy • Uflf = cmlW^j ■ n/|/ = 2(V/Ai x V/A^) • n/|/ = 

When / is a parallelogram, denote by v^, Vj and v; the vertices on / such that 


df = { eij,ejk,eki,eii}. Then by the definition of Wij and Lemma 2.3 
(curl W„) - 11 / 1 /= curl(Wij + • n/|/ 

= curl ^ 2 ^^" + ^ik + Wii) + 2 + ^kj + ' ri/l/ 

- ^curl(Wii + Wfej) • n/l/ 

= {F/fXi X V/^s) ■ n/l/ + X V/Aj) • n/l/ 

se{i,j,k,l} se{i,j,k,l} 

+ (V/^z X V/Ai + V/Aj X V/Afe) • n/l/ 

= 0 + 0 +^. 

In the above we have used Wa = Wjj = 0 and X]sG{i j k z} = 0. 

For Bij G —df, one just needs to change the sign. This completes the proof of 
the lemma. □ 


Finally, on each f G P, define 


Wf= Y. 

BijGdf 
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By the definition and Lemma |4.5[ we clearly have 
(4.4) 


f ill if / is a triangle, 

{cuTlWf)-nf\f = 

Ul it / IS a parallelogram. 
Moreover, let f'GJ-he another face of T that is different from /, then 


.X . s I j-Tjn if/,/'share an edge, 

(4.5) (curl4L/)-n/-|// = <^ 1^1 , 

10 if /, / do not share edge. 

4.2. Discrete space and Basis function. Now we are able to construct spaces 
A4A^(T) and MK^{T). It is very tempting to use Wij, for all G £■+ as a set of 
basis for the iL(curl) finite element space M.h}{T). However, the biggest problem 
of doing so is that, we are not sure whether VXi G span{Wij, for G £~*~ } or not, 
and thus can not ensure the VA4A°(T) C M.h}{T) part in the sequence (1.5). 


To ensure the exactness of sequence (1.5), similar to the 2D case, we will try 


MK^iT) = VMA°(T) e'H = span{VAi, i = 1,..., n} © H, 

MA^{T) = curl-H © (div^')®, 

where % is & space orthogonal to spanjVAi, i = l,...,n}. Again, in practice, 
it is very hard to construct orthogonal basis. Thus we relax the orthogonality 
requirement a little bit and replace © by +. Similar to (3.1), we construct the 
following: 

(4.6) AIA^(T) = span{VXi, i = 1,.. . ,n} + span{Wf, f G J-}, 

(4.7) M.A^{T) = cuT\span{Wf, f G if} + span{x — x*} 

= span{cuTlWf, f G T} + spanjx — x*}, 

where x* is a chosen point inside T. Of course this is just the construction. We 
still need to show that is exact under this construction. 

By definition, we have K G VA1A°(T), VA1A°(T) C AIA^(T), curlAlA^(T) C 
A4A^(T) and divAlA^(T) = M. Moreover, it is clear that curlAlA^(r) nspanjx — 
X*} = {0}. Th ese e stablish the exactness at the AlA°(r) and the AiA'^{T) nodes. 
To show that (1.5) is exact at the A4A^(T) node, we only need to prove that 


no none-zero vector in span{Wf, f G F} is curl free. This can indeed be done 
by counting dimensions, i.e., we will prove that the dimensions of AiA^{T) and 
MA^{T) are exactly and ^F, as indicated in (1.5). These dimensions are 
computed by explicitly constructing basis functions, as shown in the following two 
lemmas. We postpone the proof of these two lemmas to Appendix [B| 

Lemma 4.6. There exists a computable basis {q/, for f G for A4A^(T) defined 

in such that on each f G F, 


q/ -n/'l/' = 


1 *// = /', 

0 otherwise. 


Therefore the dimension of AiA^{T) is equal to the number effaces ofT. 
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Lemma 4.7. There exists a computable basis {pe, for e G for A4A^{T) defined 
in 1 4-6), such that on each e' G , 


p> — 


ife = e', 
otherwise. 


Therefore the dimension of A4A^{T) is equal to the number of edges ofT. 


Remark 4.8. By the definitions of AiA^{T) and AAA^{T), lemmas PFtI and by 
counting the dimensions, we know that (1.5) is an exact sequence. 


Remark 4.9. Lemma 4.6 indicates that for all q G A4A^(r), q • n is piecewise 
constant on the surface of T. Moreover, the normal components on faces of T 
form a unisolvent set of degrees of freedom for AiA^{T), which allows one to build 
iL(div) conforming finite element space using AiA^{T). 


Remark 4.10. Similarly, Lemma 4.7 indicates that for all q G A4A^(T), q • t is 


piecewise constant on the skeleton of T, i.e., the collection of all edges in £. More¬ 
over, the tangential components on edges of T form a unisolvent set of degrees of 
freedom for AiA^(T). However, this is not enough for building iL(curl) conforming 
finite element space, as iL(curl) conforming requires the tangential components on 
all faces, not only on edges, to be continuous across elements. 


Next, we show that the basis Pe also provides tangential continuity across faces. 
For each p G M.A^{T), its value on a face f G F can be split into two orthogonal 
parts 

Pl/ =T/(p) +Af/(p), 

where 7/(p) and A//(p) are the vector projections of p|/ onto / and its normal 
direction, respectively. We also denote by 7aT(p) the patching of 7/(p) over all 
fGF. 

By the definition of MA^{T) and Wf, it is clear that 

AAA^iT) C span{A/Xi, i = 1,... ,n} -I- span{Wij, etj G £}. 

But in general, we do not know whether VAi = — is in span{Wij^ G £} 

or not. However, if only considering the tangential component, one has the following 
nice property: 


Lemma 4.11. Let Cij G and Pe.^. be the basis function of AAA^(T) associated 
with edge e^-. Then 


TdriPeij) =TdTi\eij\Wij). 


Proof. Clearly, Tf(Wij) is nonzero on / only when both and Vj lie on /. Note 
that Equation (3.3) is still true in 3D. Therefore, one has 


raT(VA,) = - Tot{W,,) 

j such that 
Cij G £ U £f 


- E raAmj). 

j such that 
Cij G £ 


In the above we have used the definition of Wij to cancel out terms on e^i G 
£f (if there exists any) that are not connected to vertex v^. Hence £ 
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span{TdT{Wij), Cij S £}, which together with the definitions of Wf and 
further implies that 

TdT{MA^{T)) C spaniTdriWij), Cij G £} 

= span{TdTiWtj), Cij e £+}. 

Therefore, by Lemma |4.3| and by comparing the tangential components on each 
edge, one must have TdriPaj) = TdT{\eij\Wij) for all S £+. This completes the 
proof of the lemma. □ 


Remark 4.12. Lemma [4.11| tells us that the tangential component of each basis 
function on dT is completely determined by the tangential component of Wij. 
Let T and T' be two polyhedra sharing a face /, and let e^- be an edge of the 
polygon /. Then, by the definition of Wij and Assumption 1, we know that 
has continuous tangential component across the face /. Thus one can build iL(curl) 
conforming finite element spaces using AiA^{T). 


Remark 4.13. If £i = 0, then similar to the proof of Lemma [4.11 one can show 
VAi S span{Wij, Cij G £} and consequently Pe^. = \eij\Wij. Examples of poly¬ 
hedra with £/ = 0 include tetrahedra, pyramids and triangular prisms, but not 
rectangular boxes. 


Next, we briefly show that A4A^(T) C WA^{T) for fc = 1,2. By Equation (3.3) 
and the d efini tion of Wf, one immediately has AiA^{T) C >VA^(T). Similarly, by 
Equation (4.1) and the definition of Wf, one gets curlAIA^(r) C WA^iT). We also 
know from Equation (1.4) that x — x* € WA^iT). Combining the above with the 
definition oiMA^iT) gives MA^{T) C WA^{T). 

Einally, to ensure the approximation property of AAA^iT), for fc = 1, 2, we would 
like to have Vi A^{T) C AiA^{T). This is not easy to prove, and so far we do not 
even know whether it is in general true or not. Fortunately, we are able to prove 
this for two special types of polyhedra: 


Type I: Polyhedra with £/ = 0; 

Type II: Polyhedra with a center Xc such that for each vertex v^, 1 < i < n, 
one has 


(4.8) (Xc - V,) X ^ Xy = 0. 

3, eij^£ 

This is equivalent to say the barycenter of the point set {vj, for all G £} 
lies in the line passing through and Xc. 

The proof of the following Lemma will be given in Appendix [C| 

Lemma 4.14. On Type I and II polyhedra, one has Vi A^{T) C AIA^(r) for 

fc= 1,2. 


Remark 4.15. Type 1 polyhedra include all tetrahedra, pyramids, and triangu¬ 
lar prisms. Type II polyhedra include all parallelepipeds, all regular n-gon based 
bipyramids, the regular octahedron, the regular icosahedron, and some Catalan 
solids. 


Remark 4.16. From Lemma 4.14[ we know that for Type 1 and II polyhedra, the 
definition of MA^{T) is independent of the choice of x*, because C Vf A^iT). 
But so far we do not know whether the definitions of AiA^{T) and AiA^{T) are 
independent of the linear combination given in Equation (14.31) or not. 
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Remark 4.17. One may alternatively define an i7(curl) conforming finite element 


M.h}(T) = span{Wij, for e^- 




which contains A^(T) according to Lemma 4.4 for all polyhedra satisfying As¬ 


sump tions 1-2 (not restricted to Type 1 and 11 polyhedra). Moreover, by Remark 
4.13 it is clear that A4A^{T) = M.A^{T) on Type 1 polyhedra. However, as men¬ 


tioned in the beginning of this section, in general we do not know whether the 


alternative construction fits into a discrete exact sequence similar to (1.5) or not. 


4.3. Examples. We show that our construction reproduces known iJ(curl) and 
i7(div) elements on tetrahedra, rectangular boxes, pyramids, and triangular prisms. 
Then, we shall construct elements on a regular octahedron, which has never been 
done before. 

In the construction, basis functions are computed according to the proof of lem¬ 
mas and which is given in Appendix]^ Wachspress coordinates are used to 
define A^. The computation can be done using any computer algebra system. The 
results are listed below: 

(1) On any tetrahedron, there exists a unique set of barycentric coordinates. 

One can indeed easily prove that AIA^(T) = yVA^(T) = A^{T) for 

A: = 0,1, 2. No computation is needed. 

(2) On a rectangular box (0, hi) x (0, / 12 ) x (0, by using the standard tensor 
product basis: 


Ai — 


A 3 — 


As — 


{hi - x){h 2 - y)ih3 - z) 

xyjhs - z) 
hih2h^ 

{{hi - x){h 2 - y)z 


A 2 — 


A 4 — ~ 


^1^2^3 


A? — 


xyz 


Afi — 


As — 


x{h2 - y){hz - z) 
/ll/l2^3 

{hi - x)y{h 3 - z) 
hih2h^ 
x{h 2 - y)z 
/11/12/13 
{hi - x)yz 
AilAl2^3 


hih2h^ ’ 

Our construction gives 

AiA^{T) = Qo,i,i X Qi,o,i X Qip.o, 

AiA^{T) = Qipfl X Qo,i,o X Qo,o,i, 

where Qij.k = span{x^y^ , 0 <*</, 0<j<J, 0<A:< K}. These 
are identical to the lowest order Nedelec element defined in [5S] . 

Through the calculation, we also notice that on a rectangular box, the 
spaces AiA^{T) and A 4 A^{T) are much smaller than the spaces WA^{T) 
and >VA^(r) constructed in [24]. For example, one can easily see that 
W12 G WA^{T) but not in MA^{T). This indicates that there do exist 
redundant components in WA^{T) and WA^{T). 

(3) On a pyramid our construction is identical to the Whitney elements con¬ 
structed by Gradinaru and Hiptmair in [2fi] . if starting from the same 
AiA^{T) as in |2^. Since Sj = 0, one can use the simplification given in 
Remark |4.13 which coincides with the construction process in [55]. Thus 
we omit the details here. 
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(4) On a triangular prism with base defined by (0,0), (1,0), (0,1) and the 
vertical limits 0 < z < 1, we use the following barycentric coordinates: 

Ai = (l-a:-j/)(l-z), A2=a;(l-z), A3 = j/(l - z), 

A4 = (1 - a: - j/)z, A5 = xz, Ae = yz. 


Our construction gives 


MA^{T) 

MK^iT) 


(ai - 032/) + (04 - aQy)z 
(02 + a^x) + (05 + aGx)z 
a-j + a^x + ogz 


Oi G K for 1 < z < 9 


aix + 02 
aiy + 03 
04Z 05 


Oi G M for 1 < i < 5 




which is identical to the lowest order elements on triangular prism con¬ 
structed by Nedelec in m- 

( 5 ) Consider a regular octahedron, with vertices vi : (0, 0,—1), V 2 : (1,0,0), 
V 3 : (0,1,0), V 4 : (—1,0, 0), V 5 :, (0,-1, 0) and ve : (0,0,1). The analytical 
form of basis functions would be to complicated to be enclosed in this 
paper, or to be analyzed directly. Here we draw the graph of two basis 
functions for in Figure]^ In Matlab, we are also able to show 

that C M.h?{T) by computing certain linear combinations of the basis 
functions on a fine enough point grid, that reproduces constant vectors 
[1,0,0]*, [0,1,0]* and [0,0,1]* on all grid points. This numerically verifies 
that C MA^(T), which agrees with the theoretical result. 




Figure 5. Two basis functions for AiA^{T) on the regular octa¬ 
hedron. The normal component of the basis function is equal to 1 
on the shaded face and 0 on all other faces. 

We end this section with a brief discussion of elements on general hexahedra. 
Similar to the 2D quadrilateral case, the lowest order Raviart-Thomas element can 
be defined on hexahedra via Piola transform associated to a trilinear isomorphism, 
but requires asymptotically parallelepiped grid [S] in order to have good approxima¬ 
tion rate. More results on the general hexahedral Nedelec-Raviart-Thomas elements 
can be found in the recent work m and references therein. In 3D, it is also possible 
for the image of the cube under a trilinear isomorphism to have non-planar faces. 
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By working on the physical hexahedra directly, we can avoid this problem com¬ 
pletely. However, a general hexahedron does not satisfy Assumption 2, and does 
not belong to either T ype I or 11. Nevertheless, a quick examination shows t hat 


AiA^{T) from Remark 4.17 is still well-defined. Similar to the proof of Lemma 


4.4 


but requiring a more subtle treatment on Cki G £f, one can still show that Aih}{T) 

contains Vi A^{T) and consequently its curl contains Therefore, A4A^{T) may 
be used to build i7(curl) conforming finite element spaces on hexahedral meshes. 
In contrast, the situation for AAA^(T) is much more complicated, as we may not 
be able to keep the normal components on faces to be constants. Hence it remains 
a topic for future research. 


Appendix A. Adjacency matrices of a convex polyhedron 

For a convex polyhedron T, we introduce a few integer-valued matrices related 
to the shape of the polyhedron. For convenience, let us temporarily index the edges 
in by e^, for 1 < j < #A, and the faces in by /fc, for 1 < /c < Such kind 
of edge and face indices are only used in this section. In other parts of the paper, 
we do not index edges or faces of a polyhedron T by a single integer, in order not 
to be confused with the integer indices for vertices. 

Define matrices 




1 

fl 

if Ci G dfj 

j^FtoE . 


such that = < 

-1 

if e* G -dfj , 



1 

[o 

otherwise 



1 

r-1 

if Ci starts from v^- 

j^VtoE . 


such that AY/°^ = < 

1 

if Ci ends at Vj 



1 

[o 

otherwise 


For each face fi £ denote by n{fi) the number of edges in fi. For each 

Vi G V, denote by n(vi) the number of edges connected to v^. Define = 
(^A^toEyj^FtoE g l^#Fx#F ^ ^^VtoEy^VtoE g K#Vx#V^ 

to see that the entries of and are 

f n{fi) if i = j, 

M£ = < — 1 if /i) /j share an edge, 

[ 0 otherwise, 

and 

{ n(v,) iff=i, 

— I if Vi, Vj are connected by an edge, 

0 otherwise. 

To study the rank of and , let us first state a well-known result: 

Lemma A.l. Let M be an irreducible and (weakly) diagonally dominant square 
matrix, then M either has full rank or a rank 1 deficiency. 

Proof. For reader’s convenience, we provide a brief proof below. A square matrix is 
called irreducibly diagonally dominant if it is irreducible, weakly diagonally domi¬ 
nant but in at least one row is strictly diagonally dominant. Irreducibly diagonally 
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dominant matrices are non-singular. Now, by changing only one entry in any cho¬ 
sen row of M, we can make it irreducibly diagonally dominant. Since changing one 
row of a matrix can at most modify its rank by 1, therefore M must either have 
full rank or a rank 1 dehciency. □ 


Then, we have 

Lemma A. 2 . Matrix has rank {4!^ F—1), and Ker{M^) = span{[l^\ ... ^1]*}. 


Proof. By using the adjacency graph of the faces of T, it is not hard to see that 
is irreducible. Since the number of faces adjacent to each given face fj is e qual 
to n{fi), we know that is weakly diagonally dominant. By Lemma A.l 
either has full rank or a rank 1 deficiency. Indeed, M has a rank 1 deficiency, since 
one can explicitly compute that [1,1..., 1]‘ G Ker{M^). This completes the proof 


of the lemma. 


□ 


Lemma A.3. Matrix has rank {ffV — 1), and Ker{M'^) = span{[l,l... ,lY}. 

Proof. The proof is similar to the proof of Lemma |A.2[ □ 

Finally, we mention another important property of the adjacency matrices: 
Lemma A.4. It holds that 

(A.l) = 0 and = 0. 

Indeed, we have 

= range{A^^°^) and Ker{{A^*°^Y) = range{A^*°^). 

Proof. By using the adjacency relations, it is elementary to prove ( A.l[ ). Conse¬ 
quently, one has 

range{A^*°^) C and range{A^*°^) C 

Now, by lemmas A. IH we have rank{A'^*°^) = ffV — 1 and rank{A^*°^) = 
ffF — 1. The lemma follows immediately from using the rank-nullity theorem and 
counting the dimensions. □ 


Appendix B. Proof of lemmas [L6l and [TTI 
To prove Lemma |4.6[ we first denote 

q/ = c/,o(x - X*) -I- ^ Cf jcxiAWp 


and then show that there exists Cj j, for / G F} such that qj satisfies Lemma 

Denote by d/ the distance from x* to face /, and by \Tf\ = \df\f\ the volume 


4.6 


of the pyramid with base / and apex x*. For convenience, denote 


5 


f.f = 


if / = f, 

otherwise. 


For each /' G F, denote by F{f') the set of all faces in F that share an edge with 
/'. Clearly, the number of faces in F{f') is equal to the number of edges of polygon 
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/', which is denote by n{f'). Then, on each /' G F, we want {c/_ 0 j Cy for / G 
to satisfy 


<^/J' = q/ • n/'l/' = c/.o(x - X*) • ^ c^_;curlfTj • nf,\f, 

n(f') Cy ? 

= c/.od/' + Cfj, -jjjj- - XI ITT’ 


/ 6 ^(/') 


I/' 


where in the last step we have used equations (4.4)-(4.5). Multiplying both sides 
of (B.l) by I/'I and sum up over all J'gF gives 

1/1= E c/,od/'|/'l+0 = 3c^,o|r|, 


which implies 


c/.o = 


3|T|' 


Now, Equation (B.l) can be rewritten into, for each /' G F, 


i(/')c/,/'- E '^/./= ^/J'l/'l- 

feFf) 


\TfL 

\T\ 


l/l- 


This provides a linear system for solving Cf f, for all f G F, where the coefficient 
matrix is exactly defined in Appendix]^ Note the right-hand side of the above 
linear system is obviously orthogonal to Ker{M^), as 

eEi/i-^i/i)=o. 

Therefore the linear system is solvable. This establishes the existence of q/ satisfy¬ 
ing qj -rifilfi = 6fji. From the construction we also know that q/ is computable, 
with details given at the end of this section. Moreover, q/ is indeed uniquely de¬ 
fined since by setting j = 1 for all f G F, i.e., by making the coefficients in 
Ker{M^), one would get 

E Cy ^curl = curl E = curlO = 0, 

where we have used the simple fact that — 0 according to the definition 

of IT/. 

It is not hard to see that {q/, for / G F} is linearly independent. Again, by 
using ~ have 

dimA4A^(T) < dim curl (^span{Wf, f G -|- dim span{x — x*} 

< dimspanjlT/, / G F} + 1 

< (#F - 1) + 1 = #E. 


Combining the above, {q/, for / G F} must form a basis for M h?{T ) and conse¬ 
quently dim A4A^(T) = ^F. This completes the proof of Lemma 


4.6 
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Next we prove Lemma [4 .71 The idea is similar to the proof of Lemma [4.6| We 
express 

n 

Pe ~ ^ ^ ^ ^ ’ 

i=l /6JP 

Now, let e' G Denote by v„ and the starting and ending vertices of e', and 
by fi/fr the faces to t he le ft/right of edge e', seeing from outside of T. Then, by 
Assumption 1, Lemma 


4.3 


and the definition of W/, one has 


^e,e' — Pe ‘ fe 


-Oe 


./ — ^ ^ ' f e^ 

ae,/3 ^e,fi ^e^f-r 


/e^ 


|e'| |e'| 

which we further rewrite into 

(ld.2) f^e,a “t“ ^e,/3 “t“ ^e,f[ ^e,/,. — ^e.e^C |- 

The above equation holds on every e' € and thus gives us a linear system with 


#£1 equations and #V + 4j^F = + 2 unknowns. Denote by A : 

the coefficient matrix of this linear system. It is not hard to see that, under proper 
ordering, one has 

where and A^*°^ are as defined in Appendix |a| 

By Lemma |A.4[ we have 


AM = 


0 


g ]j(#iJ+2)x(#i5+2)^ 


Consequently, by lemmas we know that rank(A) = rank{A*A) = (#D — 

1) + (#i^ — 1) = #E and Ker{A) is spanned by the following two vectors: 


(B.3) 


[1,1,..., 1,0,0,..., 0]*, with (#D) I's and (#T') O's, 
[0, 0,..., 0,1,1,..., 1]*, with (#D) O's and (#J^) I's. 


Then, the linear system (B.2) is solvable. Moreover, we realize that Pe is indeed 
uniquely defined, as all coefficients in Ker{A) only generate zero functions because 
Er=iVA. = OandE/e^^^/ = 0. 

We can similarly show that {pe, for e G £“*■} is linearly independent, and thus 
by counting dimensions, it form a basis for AtA^(T). This completes the proof of 
Lemma 14.71 

Finally, we briefly discuss how to compute the basis functions in practice. Using 
elementary linear algebra, it is not hard to see that: 

(1) To compute q/, one needs to solve a linear system M^u = b, where G 
has a non-trivial kernel containing all constant vectors, and b G 
Ker{M^)^ = Range{M^). Indeed, solving M^u = b is equivalent to 
solving a non-singular square system 


'M^ 

1 


U 


b 

1‘ 

0 


0 


0 


where 1 denote a constant column vector with all entries equal to 1. 
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(2) To compute qe, one needs to solve a linear system Au = b, where A = 
^j^vtoE ^FtoEj g ]j#Ex(#E+ 2 ) kernel spanned by vectors 

in ( |B.3[ ). Indeed, solving Alu = b is equivalent to solving a non-singular 
square system 


~ j^VtoE 

j^FtoE' 


h 

1* 

0‘ 

U = 

0 

0* 

i‘ 


0 


Appendix C. Proof of Lemma [4.141 
For Type I polyhedra, the proof is easy. By Remark 


on each S . Thus by Lemma 


4.4 


4.13 


we have pe,^- = \eij\Wij 
one immediately gets Pf A^(r) C AIA^(T). 


This, together with the fact that curl(a x x) = 2a for all a S implies that C 
AIA^(T). Finally, since span{x — x*} C AIA^(T), we have Pf A^(T) C A4A^(T). 
Now let us consider Type II polyhedra. From the proof of Lemma |4.4[ one has 

((a X (vi - Xc)) • ^ ((a X (v^ - x^)) • Tij)Wij = 2a x (x - x^). 


:i,G£ + 


Cij^S 


for all a S 

(C.l) 


If we can show that 
((a X (vj - Xc)) ■ 


eij&E + 


^ CfWf G MA\T), 

/GF 


this together with the face that C VAIA°(T) C AIA^(T) will imply Vi A^{T) C 
A4A^(T). And consequently one will be able to prove that V^A'^iT) C M.A^{T). 
Next, we focus on prove the existence of a set of coefficients {C/, for / G F} that 
satisfies ( |C.1[ ). 

For each e^- G , denote by and the faces on the left and right side of 
respe ctively, seeing from outside of T. Notice that the right-hand side of Equation 
(C.l) can further be written into ~ Thus it remains to 


prove that the system 

(C.2) Cfij - Cfij = (a X (vi - Xc)) 


for all Cij G 


is solvable. The coefficient matrix of system (C.21 is exactly , as defined in 

Appendix]^ Denote the right-hand side vector of system (C.2) by 

b = [(a X (v, - Xc)) • (,£+ G 


The linear system (C.2) is solvable only if 

b G Range{A^^°^) = Ker{{A^*°^ f)^. 

By Lemma [AaI we have Ker{{A^^°’^Y)^ = range{A^*°^)^ = Ker{{A^^°^Y). 


Therefore, System (C.2) is solvable as long as (A^*°^)*b = 0 , which can be explic¬ 
itly written as 

- E 

j,ejie£ + 

where we conveniently denote by bij the entry of vector b corresponding to e^- G £'^ . 


(C.3) 


E ' 

i,eoG£ + 


bji — 0, 


for all 1 < i < 


According to (C.21, bij can be viewed as the jump of coefficient C/ across the edge 


Cij. Thus the constraints given by (C.3) are equivalent to say that, the summation 


of such jumps over all edges connecting to one given vertex should be 0. By the 
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definition of b and the fact that x Tij = vj x Equation (C.3) is equivalent 
to 

E (ax(vi-Xc))-Xy- ^ (ax(vj-Xc))-Tji = ^ (ax (v,-Xc))-ry = 0, 

3,eij^£+ j,ejiG£+ 

for all 1 < i < n, which is true on Type II polyhedra. In other words, we have 
shown that for Type II polyhedra, Equation (C.2l is solvable. This completes the 
proof of Lemma [4.14| 
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